/*
 * =====================================================================================
 *
 *       Filename:  process.h
 *
 *    Description:  matrix process (like filter)
 *
 *        Version:  1.0
 *        Created:  2009年06月07日 09时49分18秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang\nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */
#ifndef PROCESS_H
#define PROCESS_H

#include "ncvector.h"
#include "ncmatrix.h"

/**
 * \brief Image filter with the kernel template h
 *
 * \param mat is the input data which should be filtered
 * \param h is the kernel template
 *
 * \retrun return the data after filter.
 */
template<typename T>
inline NCmatrix<T> filter(const NCmatrix<T> &mat, const NCmatrix<T> &h)
{
	int i,j, hi, hj, ii,jj,row = mat.row(), column = mat.column(), hrow=h.row(), hcolumn=h.column();
	int kcr = (hrow-1)/2, kcc=(hcolumn-1)/2;
	NCmatrix<T> result(row,column);
	double tmp;
	for ( i=0; i<row; i++ )
	{
		for ( j=0; j<column; j++ )
		{
			tmp = 0.0;
			for ( hi=0; hi<hrow; hi++ )
			{
				for( hj=0; hj<hcolumn; hj++ )
				{
					ii = i-kcr+hi;
					jj = j-kcc+hj;
					if( ii>=0 && ii<row && jj>=0 && jj<column )
					{
						tmp += h[hi][hj] * mat[i-(hrow-1)/2+hi][j-(hcolumn-1)/2+hj];
					}
				}
			}
			result[i][j] = tmp;	
		}
	}
	return result;
}

/**
 * \brief  Generate gauss kernel template
 * 
 * \param row
 * \param column
 * \param sigma
 *
 * \return return the gauss kernel template	
 */
inline NCmatrix<double> gauss (int row, int column, double sigma)
{
	NCmatrix<double> result(row,column);
	int i,j;
	double tmp=0,crow=(row-1)/2.0, ccolumn = (column-1)/2.0;
	for( i=0; i<row; i++ )
	{
		for( j=0; j<column; j++ )
		{
			result[i][j] = exp( - ((i-crow)*(i-crow) + (j-ccolumn)*(j-ccolumn))/(2*sigma*sigma) );
			tmp += result[i][j];
		}
	}
	for ( i=0; i<row; i++ )
		for ( j=0; j<column; j++ )
			result[i][j] /= tmp;
	return result;
}

/**
 * \brief Generate difference of gauss kernel template D(x,y,sigma) = G(x,y,k*sigma)-G(x,y,sigma) 
 *
 *
 */
NCmatrix<double> dog (int row, int column, double sigma1, double sigma2)
{
	NCmatrix<double> result(row,column), g1(row,column), g2(row,column);
	int i,j;
	double tmp1=0,tmp2=0,crow = (row-1)/2.0, ccolumn = (column-1)/2.0;
	double dis;
	for(i=0;i<row;i++)
	{
		for(j=0;j<column;j++)
		{
			dis = (i-crow)*(i-crow) + (j-ccolumn)*(j-ccolumn);
			g1[i][j] = exp(-dis/(2*sigma1*sigma1));
			g2[i][j] = exp(-dis/(2*sigma2*sigma2));
			tmp1 += g1[i][j];
			tmp2 += g2[i][j];
		}
	}
	for(i=0;i<row;i++)
	{
		for(j=0;j<column;j++)
		{
			result[i][j] = g1[i][j]/tmp1 - g2[i][j]/tmp2;
		}
	}

	return result;
}
/**
 * \brief laplacian of Gaussian
 *
 *
 */
NCmatrix<double> log (int row, int column, double sigma)
{
	NCmatrix<double> result(row, column);
	int i,j;
	double crow = (row-1)/2.0, ccolumn = (column-1)/2,sigma2=sigma*sigma,sigma4=sigma2*sigma2;
	double dis;
	for ( i=0; i<row; i++ )
	{
		for( j=0; j<column; j++)
		{
			dis = (i-crow)*(i-crow) + (j-ccolumn) * (j-ccolumn);
			result[i][j] = (dis/sigma4-2/sigma2)*exp(-dis/sigma2);
		}
	}
	return result;

}

/**
 * \brief x-direction gradient of matrix X
 *
 */
template<typename T>
NCmatrix<T> gradientx(const NCmatrix<T> &X)
{
	int i,j,row=X.row(), column=X.column();
	NCmatrix<T> result(row,column);
	
	for( j=1; j<column-1; j++ )
	{
		for( i=0; i<row; i++ )
		{
			result[i][j] = (X[i][j+1]-X[i][j-1])/2.0;
		}
	}

	for( i=0; i<row; i++ )
	{
		result[i][0] = X[i][1]-X[i][0];
		result[i][column-1] = X[i][column-1]-X[i][column-2];
	}
	return result;
}

/**
 * \brief y-direction gradient of matrix X
 *
 */
template<typename T>
NCmatrix<T> gradienty(const NCmatrix<T> &X)
{
	int i,j, row=X.row(), column=X.column();
	NCmatrix<T> result(row,column);
	for(i=1;i<row-1;i++)
	{
		for(j=0;j<column;j++)
		{
			result[i][j] =  (X[i+1][j]-X[i-1][j])/2.0;
		}
	}
	for( j=0; j<column; j++ )
	{
		result[0][j] = X[1][j] - X[0][j];
		result[row-1][j] = X[row-1][j] - X[row-2][j];
	}
	return result;
}
#endif
